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The problem of a physical relevance (meaning) of percolation in supercritical fluids is addressed considering a 
primitive model of water. Two different criteria, physical and configurational, are used for the cluster definition 
in Monte Carlo simulations over a range of pressures to determine the percolation line and skewness, and a 
theoretical analytic equation of state is used to evaluate response functions. It is found that both criteria yield 
practically the same percolation line. However, unlike the findings for simple fluids, the loci of the response 
function extrema exhibit density/pressure dependence quite different from that of the percolation line. The 
only potential coincidence between the loci of the extrema of a thermodynamic property and a detectable 
structural change is found for the coefficient of isothermal compressibility and Voronoi neighbors distribution 
skewness maximum. 
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1. Introduction 

It has been well known that at certain thermodynamic conditions some properties of fluids may ex- 
hibit a rather sudden sharp change. Typical examples are a rapid change of electric conductivity of aque- 
ous solutions |1, 2], metal-nonmetal transition in supercritical fluids I3J, or a sudden change of electric 
and dielectric properties of composites |4]. Formation of clusters (aggregates) of molecules and/or a net- 
work has been made responsible for these observations. 

Unlike the crystalline phase, studies of percolation in a fluid phase encounter a number of both tech- 
nical and fundamental (methodological) problems. First of all, in continuous systems with a continu- 
ous interaction potential it is not clear how to define a bond. Hence, several criteria have been put for- 
ward 1 5, 6]. Similarly, using molecular simulations on finite samples the question arises how to define and 
detect the presence of an infinite cluster. Consequently, practically no theoretical results are available for 
fluids while various conjectures on percolation in fluids represent a mere extension of the results orig- 
inally obtained for lattice systems. We have addressed these problems in a series of recent papers (7|-|9|] 
with the following two most important results. It was shown that (i) unlike the lattice systems, the percola- 
tion in continuous systems does not follow any scaling law and cannot be thus characterized by universal 
constants and that (ii) the wrapping probability (the cluster wraps the system if, starting from any parti- 
cle of the cluster and moving along interparticle bonds, there is a possibility of getting to an image of that 
particle in another replica 1 10]) provides a more accurate estimate of the percolation threshold compared 
to the crossing probability (the cluster crosses the system if the maximal distance between some pairs of 
its particles is equal to or greater than the box length) 

Instances of the observed sharp changes in some properties of fluids usually refer to mixtures. How- 
ever, it is also well known that the response functions (derivatives of thermodynamic potentials) of pure 
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fluids exhibit extrema in the supercritical region fllL It means that some thermodynamic properties 
change the course of their dependence on the thermodynamic conditions and the natural question to be 
asked is what is behind this change in the behavior at the molecular level. A potential answer is that these 
changes may be associated with the formation of an infinite cluster. If this is the case, then the percolation 
line, i.e., a line connecting all state points where the percolation transition/threshold occurs, becomes a 
quantity with a real physical content and will not represent a mere mathematical construction extending 
otherwise meaningful results. 

In a recent paper 1 12] we have addressed this important issue by studying (by Monte Carlo simu- 
lations) the percolation and various structural characteristics of two simple fluids, the square-well and 
Lennard-Jones ones. We have found that out of all the considered extrema lines of response functions 
(RFEL; i.e., the lines connecting the state points at which the given response function, as a function of 
temperature T at constant pressure P, reaches its extremum) only the RFEL of the coefficient of isother- 
mal compressibility has closely followed the course of the skewness extremum line (SEL; i.e., the line 
connecting the state points at which the skewness of distribution of the number of nearest neighbors 
in terms of the Voronoi tessellation flallal reaches its maximum) whereas no definite conclusion has 
been possible to draw concerning other response functions. A similar study has been conducted recently 
by Idrissi et al. 1 16] who investigated the relation of thermodynamic response functions with the Voronoi 
tessellation properties in supercritical ammonia. They found a potential relation between the thermal ex- 
pansion coefficient and other Voronoi tessellation properties, e.g., the volume distributions. Specifically, 
the volume distributions broaden near the temperature where the thermal expansion coefficient reaches 
its maximum. Nonetheless, the relation between the percolation transition on one hand and the thermo- 
dynamic response functions and/or Voronoi tessellation on the other hand has not been fully explored 
yet. 

Accounting for a number of uncertainties when dealing with simple fluids the inconclusive results 
mentioned above are not surprising. As a more severe test of the original surmise of a coincidence of 
the percolation and SEL lines and the RFEL lines, a system with a natural formation of chains or clusters 
should be considered. Water is a typical example of a system with association and we have therefore 
carried out the same project considering a simple model of water, i.e., the primitive model, for which (i) 
the criterion of bonding is unique and (ii) for which an analytic equation of state is also available Il7ll . 

In the present paper we report the results for the same properties as in the previous study 1 12]. In the 
next section we provide all necessary definitions and computational details and in section 3 we present 
and discuss the obtained results. 



2. Basic definitions and computational details 

In this paper we consider a primitive model of water (PMW) descending from the realistic parent 
TIP4P model hereafter referred to as the PM/TIP4P model 1 18]. This is a simple short-range archetype of 
real water faithfully reproducing its structure and results from a non-speculative modeling applied to 
the parent model 1 19] . In the approach used, the Lennard-Jones interaction in the realistic parent model 
is replaced by a hard-sphere repulsion, while the repulsive and attractive electrostatic interactions at 
short separations are approximated by hard sphere repulsion and square-well attraction, respectively 
(for details see (l^|20|]). Thus, the functional form of the model reads as 

w PM (l,2) = MHs(roo.doo) + £ umirij-,dij)+ £ wswtoy) 

ie{l},;e{2} ie{l},i'e{2} 
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and the summation in the second term of l[T) runs over the pairs of the like sites and in the third term of 
the unlike sites. The above PMW was developed to serve as a short-range reference system for developing 
a full perturbation theory of water (illl . It was shown to faithfully reproduce the structure of water as 
well as to capture a good deal of its thermodynamic properties including the observed anomalies. The 
parameters of the model are presented in table [l] It is worth recalling that the potential «phb in (1) 
incorporates all repulsive interactions of the PMW; after switching off all attractive interactions in the 
PMW we are left with a rather strange hard body called pseudo-hard body (PHB) (22)]. It is not a simple 
hard body because it possesses a flavor of nonadditivity and currently no theory for the fluid of PHB's is 
available lH. 



Table 1. Parameters of the considered PM/TIP4P model of water 1 18]. 
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For the PMW with its strongly orientation dependent short-range attractive interaction, it is only nat- 
ural to consider two molecules to be bonded when their interaction energy uij is -£"hb (the parameters 
of the model satisfy the so-called conditions of saturation which means that no more than one bond 
can be established between a pair of molecules). The other criterion used to define a bond between two 
molecules is the energy-based Hill's definition according to which two molecules are bonded when their 
relative kinetic energy is less than their interaction energy (HHil . Specifically, if v; and vj are the velocity 
vectors of particles i and /, respectively, and m is their mass, then they are considered to be bonded if 

-(vy-Vj) ^- Uij . (4) 

We used the common Metropolis Monte Carlo (MC) simulation in an NVT ensemble (zsfa both for the 
location of the percolation line and SEL which means that two types of simulations had to be performed. 
The wrapping probability R w as a function of reduced number density p* at a given temperature T 
and number of particles N have to be determined first. Then, the percolation threshold density p* c {T) 
is calculated as the density where R w at given T for different N intersect in one point (more accurate 
numerical procedure is described in |9]). Details of this simulation were as follows: (i) Wwas set to 1000, 
2000, or 4000, (ii) each R w {p*,{N, T\ - const) was determined for at least 22 (with maximum up to 42) 
densities p* with step 0.002 such that the minimum of the mean of that R w was less than 0.03 and its 
maximum was greater than 0.97, (iii) the systems were equilibrated by performing at least 32 x 10 4 AT 
MC steps, (iv) equilibrium configurations used for the detection of the presence of the wrapping cluster 
were separated by 37V MC steps, (v) at least 7 x 10 6 configurations (up to 124 x 10 6 ; increasing with T 
going closer to critical point and with smaller N) were investigated to reach a sufficiently low error of 
measurement, and (vi) to apply the energy criterion, the velocity components were randomly assigned to 
every particle from the Gaussian distribution characterized by the temperature of the system (for further 
details see the previous paper |9]). The simulation details for determining the SEL were as follows: (i) N 
was set to 1000, (ii) each skewness j\{p* ,T - const) was determined for at least 28 values of p* (with step 
0.002), (iii) the systems were equilibrated by performing at least 10 5 N MC steps, and (iv) 10 6 equilibrium 
configurations, separated by 100AT MC steps, were used for the Voronoi polyhedra analysis in which the 
water molecules were represented only by the position of their oxygen-like sites. The parameters of all 
the simulations mentioned above were set so as to maintain the acceptance ratio around 1/3. 

The response functions that we consider to be of interest are the coefficient of thermal expansion, 
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the coefficient of isothermal compressibility, 

V \dPjj p \dP) T 

and the isobaric heat capacity, 

1 (dH\ 

where H is enthalpy Finally, an extension of the vapor-liquid (V-L) coexistence curve to the supercritical 
temperatures was determined by extrapolating the Clausius-Clapeyron equation \r\P - a-bIT, where a 
and b are fitted parameters for the model considered. Both the response functions (and corresponding 
RFEL's) and the extension of V-L line were determined numerically from Vlcek and Nezbeda equation 
of state {vh . An analytic theory of the considered PMW was developed using a variant of the 2nd or- 
der thermodynamic perturbation theory (TPT2) 12611 which was shown to be very accurate. In general, 
the residual Helmholtz energy per particle, a m& - A res /N for the fluids consisting of molecules with m 
hydrogen-like sites and n complementary iV-sites assumes the form jl7ll2lll 

pXm = j8a^ B + m(l - v) +ln x , (8) 

where xq is given, using the reduced density, p* = P^q , by 

v m 

x o = 7 ; 77T - O) 

[1 + mp*vl+- + (mp* vyi + -+\ 

and parameter v, measuring 'non-saturation', satisfies the following cubic equation: 

m 2 p* 2 I+-+v 3 + {{2mn-m 2 )p* 2 I + -+ + mp*I + -jv 2 

+ [(«-m)p*7 + _ + l]v-l = 0. (10) 

Quantities I are integrals involving the correlation function of the underlying PHB fluid and are avail- 
able in a parametrized form (for further details and values of the parameters see reference 1 21]). For 
compressibility factor, z = /3PI p, one gets the following result 

= zpHB + p^[m(l-v)+lnx ]j . (11) 



3. Results and discussion 



Following the strategy outlined above using simulations we determined two types of the percola- 
tion line corresponding to the two different definitions of bonds (configurational and physical) and the 
results are shown in figures [l] and [2] The first finding is discernible. The configurational and physical 
criteria yield the identical percolation line. In our previous papers we always used these two criteria be- 
cause they always produced different lines giving rise to an ambiguity in the interpretation of the results. 
The reason why this is not the case for the considered model probably lies in the choice of thermody- 
namic conditions: unlike in the previous studies, we operate in the low temperature supercritical region. 
Stillinger's criterion imposes that utj < 0. Hill's criterion is stronger, it requires that Uij + ktj < where 
kij > is the relative kinetic energy. Since the two percolation lines are practically identical, it means that 
the relative kinetic energy of two molecules is negligible in comparison with the H-bonding energy. This 
argument is supported by the results obtained previously for the same model at higher temperatures. 
Due to a high temperature, the kinetic energy attains large values and starts competing with the H-bond 
energy leading to different conditions for percolation and thus to two different percolation lines. 

Unlike in our previous papers on percolation, where we investigated this phenomenon over a wide 
range of thermodynamic conditions, in this paper we have focussed on supercritical temperatures as 
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close to the critical point as possible. The still open question whether the percolation line originates/termi- 
nates in the critical point is one of other important issues. Examination of the figures leads to the result 
similar to that obtained for simple fluids: in the P-T plane, the percolation line tends to the critical point 
as predicted by theory but this is not the case in the T — p plane. However, when discussing this issue one 
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Figure 2. The same as figure[T]in the temperature-density plane. 



23301-5 



J. Skvor, I. Nezbeda 



has to bear in mind that while the model and the range of the considered temperatures have removed 
the ambiguity concerning the definition of bonding of two molecules, a certain uncertainty yet remains. 
Namely, evaluation of the model properties from an equation of state. The used equation was shown 
to be quite accurate at high temperatures but its accuracy deteriorates with decreasing temperature. 
Particularly, the location of the critical point determined from this equation is only approximate. 

As usual, we show the percolation lines and the RFEL' in two different frames: in the P - T plane 
which is common when discussing any transition and in the T — p plane which is a common practice 
when dealing with percolation. As we see, results in both cases are a bit different and we attribute this 
difference to the use of only an approximate equation of state as mentioned above. To discuss the results 
for the response functions, it is useful to distinguish two types of the response functions: (i) one, reflecting 
the geometrical arrangement of the molecules, i.e., the coefficient of isothermal compressibility, and (ii) 
the remaining functions that involve energy. In full agreement with the previously obtained results, we 
find that the coefficient of isothermal compressibility follows the course of the SEL. This is correct if we 
examine this behavior in the P—T plane, whereas in the T—p plane there is only a rough similarity in the 
course of these two functions. As regards other RFEL's, the result is similar to that found for the simple 
fluids: the percolation line follows the course which differs from that of the RFEL, although qualitatively 
there may be found the same functional dependence. 

The final remark on the obtained results concerns an 'imaginary' extension of the vapor-liquid (V-L) 
coexistence curve also shown in the figures. The relation between the percolation line and the extended 
coexistence curve in supercritical water was discussed already by Partay and co-workers (271 |28fl. In the 
first paper they pointed out a possible coincidence of both lines and in the second one they modified their 
claim stating only that the percolation line goes through higher pressures than the extension of the V-L 
line. However, quite an opposite conclusion can be drawn from the present paper. Taking into account 
that we used a different model (primitive contrary to realistic one) and that the extended coexistence 
curve has no physical sense, such disagreement is not that surprising. Moreover, the extension of the V-L 
line was already shown 1 12] to have no relation either to any of the RFEL's or to the percolation line in 
the case of simple fluids and the same holds true for the considered model of water. 

4. Conclusions 

Water molecules naturally associate and the definition of the bond is unique for the model consid- 
ered. We have, therefore, embarked on this project with the hope to find at least some definite answers 
to the questions raised in the Introduction. It means, to find a link between the percolation and/or cer- 
tain structural characteristics of the fluid and the observed changes in the thermodynamic behavior of 
supercritical fluids. The only positive, and probably generally valid result is that the rearrangement of 
molecules under pressure and witnessed by the Voronoi neighbors distribution skewness maximum gives 
rise to an extremum of the coefficient of isothermal compressibility. Concerning other response func- 
tions, we have to admit that the achieved results turned out to be disappointing leaving the question of 
the physical sense of the percolation in supercritical fluids unanswered, the same applying to the origin 
of the observed extrema of the response functions. 
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/lima nepKOsiflLjii', <t>yHi<i^iY Biflryi<y i aHa/ii3 6araTorpaHHi/iKa 

BopOHOrO B HaflKpMTMHHiM BOfli 

fl. WkbofR I. He36efleP 

(PaKyflbTeT npwpoflHW4wx HayK, YHiBepcwTeT fl.E. nypKiHtfe, Ydi Hafl/la6eM, HecbKa Pecny6/iiKa 
iHCTHTyT (JiyHflaMeHTa/ibHi/ix ochob xiMinHnx npoL4eciB, AKafleMia HayK, npara, HecbKa Pecny6/iiKa 

npo6/ieMa c(>i3H4H0T Ba>K/inBOCTi (cyTi) nepKO/ismiTy HaflKpnTUHHUx n/iwHax floaiiflxyeTbCfl Ha npwK/iafli npuMi- 

TMBHOT MOfle/li BOflH. fl,Ba pi3HWX KpHTepi'l', cf>i3HHHWi7i i KOHCflirypaL^iMHwR, BHKOpi/lCTOByK)TbCFl fl/lfl 03Ha4eHHfl 

K/iacrrepa b MOfle/iiOBaHHi MeTOflOM MoHTe Kap/io npn pi3Hi/ix TWCKax, 3 MeTOio cnpuMaHHH nepK0/isqi4iRH0T yiiHi'f 
Ta i"i nepeKOcy. TeopeTUHHe aHa/ii™4He piBHAHHfl CTaHy b n ko p hcto By€Tbc?i fl/ifl ol^hkh cj^yHKL^M BiflryKy. 3Ha- 
fifleHO, mo o6nflBa KpnTepiT flaKDTb npaKTUHHO OflHaKOBy /imiK> nepKO/iflLiii. npoTe, Ha BiflMiHy B\p, pe3y/ibTa"riB 
Pflft npocTux n/iwHiB, reoMeTpi/i4He u\cufi tohok eKCTpeMyMiB ct>yHKL(iT BiflryKy BKa3ye Ha i^i/iKOM mwy 3a/ie- 
XHiCTb Bifl TUCKy/rycTHHH, Hi>K nmm nepKO/iflL4ii. £fli/iHe noTeHijMHe cniBnafliHHfl Mix reoMeipuHHWM Mici4eM 
tohok eKCTpeMyMiB TepMOflWHaMHHoT xapaKTepwcTHKw i cnocTepexyBaHO'i CTpyKTypHOi 3MiHn e 3HatifleHO fl/ia 
Koetfiii^ieHTa i30TepMiHHOi CTnainBOCTi i MaKCHMyMy acuMeipiT po3nofli/iy BopoHoro. 

K/iiOHOBi c/ioBa: n'lH'm nepKonnu,n, ipyHKi^i) BipryKy, jiihu BipflMa, HapKpvmAHHa Bopa, M03aiKa BopoHoro 



23301-8 



